from brian import *
from brian.library.random_processes import *
from brian.library.synapses import  *
import mydataIC as mydatic

md = mydatic.ICVars()
md.__init__()
#defaultclock.dt = 0.0025*ms

def i_inj(t):
    if t%(500*ms) > 250*ms:
        #return 1.345 *nA
        return 0.4*nA
    else:
        return 0*nA
    

locDUR =1000
#BINSIZE = 50
# Parameters 
#L=20*umetre
#dia=20*umetre
#area=dia*L*math.pi
#cpf=(1*ufarad*cm**-2)*area
#gl=(0.1*msiemens*cm**-2)*area
#glint=(0.51*msiemens*cm**-2)*area
#El=-64*mV
#Elint = -64*mV
#Ek=-80*mV
#Ena=55*mV
#gna=(30*msiemens*cm**-2)*area
#gnaint=(120*msiemens*cm**-2)*area
#gk=(36*msiemens*cm**-2)*area
#gkint=(100*msiemens*cm**-2)*area
#gnap=(0.25*msiemens*cm**-2)*area
#taumna = 0.000001*msecond
#tau_m_nap = 0.01*msecond
#vinit=-63
#taue=5*ms
#taui=5*ms
#invtau=1/taue
#Ee=-10*mV
#Ei=-70*mV
#gsynapse=50*uS*cm**-2*area
#we=0.5*nS # excitatory synaptic weight (voltage)
#wi=67*nS # inhibitory synaptic weight

def spToList(a):
    l = list()
    
    for x in a.spiketimes:
        l+=((a[x]/20*1000/md.BINSIZE)/ms).tolist()
    return l

def napvtrap(x,y):
    a=x/y
    if abs(a)<0.0001:
        return  1/(2+a)
    else:
        return 1/(exp(a)+1)
    #return 1/(exp(a)+1)
    
def napefun(x,y):
    a = x/y 
    if abs(a)<0.0001:
        return 1+a
    else:
        return exp(a)
    #return exp(a)
            
 
def m_inf_na(V):
    return napvtrap(-(V+35),7.8)   
    
def h_inf_na(V):
    return napvtrap((V+55),7)

def tau_h_na(V):
    return msecond*(30/(napefun((V+50),15)+napefun(-(V+50),16)))

def m_inf_k(V):
    return napvtrap(-(V+28),15)

def tau_m_k(V):
    return msecond*(7/(napefun(V+40,40)+napefun(-(V+40),50)))

def m_inf_nap(V):
    return napvtrap(-V-47.1,3.1)

def h_inf_nap(V):
    return napvtrap(V+59,8)

def tau_h_nap(V):
    return msecond*(2*1200/(napefun(V+59,16)+napefun(-V-59,16)))

def current_inj(t):
    if t > 0*msecond:
        return -0.5*namp
    else:
        return 0*mamp

# The model

#eqsPF+=OrnsteinUhlenbeck('I',mu=0.005*nA,sigma=0.002*nA,tau=10*ms)
eqsIN=Equations('''
dvm/dt = (i_inj(t)-(md.gna*mna*mna*mna*hna*(vm-md.Ena))-(md.gk*mk*mk*mk*mk*(vm-md.Ek))-md.gl*(vm-md.El)+ge_current)/md.cpf : mvolt
dmna/dt = (m_inf_na(vm/mvolt)-mna)/md.taumna : 1
dhna/dt = (h_inf_na(vm/mvolt)-hna)/tau_h_na(vm/mvolt) : 1
dmk/dt = (m_inf_k(vm/mvolt)-mk)/tau_m_k(vm/mvolt) : 1
''')
eqsIN+=alpha_conductance(input='ge',E=md.Ee, tau=md.taue,conductance_name='epsp') 


#FREQ=30*Hz
#spiketimes = [(0,5*ms),(1,5*ms),(2,5*ms),(3,5*ms),(4,5*ms),(5,5*ms),(6,5*ms),(7,5*ms),(8,5*ms),(9,5*ms),(10,5*ms),(11,5*ms),(12,5*ms),(13,5*ms),(14,5*ms),(15,5*ms),(16,5*ms),(17,5*ms),(18,5*ms),(19,5*ms),]

    
#SG = SpikeGeneratorGroup(20,spiketimes,period=1/FREQ)

IN=NeuronGroup(1,model=eqsIN,
threshold=EmpiricalThreshold(threshold=-20*mV),
implicit=True)


# Initialization



vin =-64
IN.vm=vin*mV
IN.mna = m_inf_na(vin)
IN.hna = h_inf_na(vin)
IN.mk = m_inf_k(vin);

# Record a few trace
traceIN=StateMonitor(IN,'vm',record=[0])

#spikesIN= SpikeMonitor(IN)

run(locDUR*msecond)
figure(1)
#plot(tracePFgi.times/ms,tracePFgi[0]/nA)
#plot(tracePFge.times/ms,tracePFge[0]/nA)
#subplot(411)
plot(traceIN.times/ms,traceIN[0]/mV)
#subplot(412)
#x = spToList(spikesIN)
#hist(x,floor(md.DUR/md.BINSIZE),range=[0,md.DUR],histtype='step')
show()